/**     ----------------------------------------------------------------
*
*   @brief this file contains extra routines needed for the main test program for sgp4.
*   @author david vallado
*
*
*    this file contains extra routines needed for the main test program for sgp4.
*    these routines are derived from the astro libraries.
*
*    companion code for fundamentals of astrodynamics and applications
*                                    2007
*                              by david vallado
*
*       (w) 719-573-2600, email dvallado@agi.com
*
*    current :
*               7 may 08  david vallado
*                           fix sgn
*    changes :
*               2 apr 07  david vallado
*                           fix jday floor and str lengths
*                           updates for constants
*              14 aug 06  david vallado
*                           original baseline
*       ----------------------------------------------------------------      */

#include "sgp4_date.h"


/** -----------------------------------------------------------------------------
*    @brief   this procedure finds the julian date given the year, month, day, and time.
*    @param year        - year                           1900 .. 2100
*    @param mon         - month                          1 .. 12
*    @param day         - day                            1 .. 28,29,30,31
*    @param hr          - universal time hour            0 .. 23
*    @param min         - universal time min             0 .. 59
*    @param sec         - universal time sec             0.0 .. 59.999
*    @param jd         output - julian date                    days from 4713 bc
*    @author david vallado
*
*                           procedure jday
*
*  this procedure finds the julian date given the year, month, day, and time.
*    the julian date is defined by each elapsed day since noon, jan 1, 4713 bc.
*
*  algorithm     : calculate the answer in one step for efficiency
*
*  author        : david vallado                  719-573-2600    1 mar 2001
*
*  inputs          description                    range / units
*    year        - year                           1900 .. 2100
*    mon         - month                          1 .. 12
*    day         - day                            1 .. 28,29,30,31
*    hr          - universal time hour            0 .. 23
*    min         - universal time min             0 .. 59
*    sec         - universal time sec             0.0 .. 59.999
*
*  outputs       :
*    jd          - julian date                    days from 4713 bc
*
*  locals        :
*    none.
*
*  coupling      :
*    none.
*
*  references    :
*    vallado       2007, 189, alg 14, ex 3-14
*
* --------------------------------------------------------------------------- */

/**
 * @brief Esta función devuelve el día juliano dado por una fecha desde el año hasta el segundo
 * @param year Año. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param mon Mes.
 * @param day Días
 * @param hr Horas
 * @param minute Minutos
 * @param sec Segundos
 * @param jd Día juliano. Perderá el valor que tuviera antes de llamar a ésta función.
 */
void    sgp4_date::jday
        (
          int year, int mon, int day, int hr, int minute, double sec,
          double& jd
        )
   {
     jd = 367.0 * year -
          floor((7 * (year + floor((mon + 9) / 12.0))) * 0.25) +
          floor( 275 * mon / 9.0 ) +
          day + 1721013.5 +
          ((sec / 60.0 + minute) / 60.0 + hr) / 24.0;  // ut in days
          // - 0.5*sgn(100.0*year + mon - 190002.5) + 0.5;
   }  // end jday

/** -----------------------------------------------------------------------------
*  @brief   this procedure converts the day of the year, days, to the equivalent month
*    day, hour, minute and second.
*    @param  year        - year                           1900 .. 2100
*    @param  days        - julian day of the year         0.0  .. 366.0
*    @param mon         output- month                          1 .. 12
*    @param  day         output- day                            1 .. 28,29,30,31
*    @param hr          output- hour                           0 .. 23
*    @param min         output- minute                         0 .. 59
*    @param sec         output- second                         0.0 .. 59.999
*    @author david vallado
*
*
*                           procedure days2mdhms
*
*  this procedure converts the day of the year, days, to the equivalent month
*    day, hour, minute and second.
*
*  algorithm     : set up array for the number of days per month
*                  find leap year - use 1900 because 2000 is a leap year
*                  loop through a temp value while the value is < the days
*                  perform int conversions to the correct day and month
*                  convert remainder into h m s using type conversions
*
*  author        : david vallado                  719-573-2600    1 mar 2001
*
*  inputs          description                    range / units
*    year        - year                           1900 .. 2100
*    days        - julian day of the year         0.0  .. 366.0
*
*  outputs       :
*    mon         - month                          1 .. 12
*    day         - day                            1 .. 28,29,30,31
*    hr          - hour                           0 .. 23
*    min         - minute                         0 .. 59
*    sec         - second                         0.0 .. 59.999
*
*  locals        :
*    dayofyr     - day of year
*    temp        - temporary extended values
*    inttemp     - temporary int value
*    i           - index
*    lmonth[12]  - int array containing the number of days per month
*
*  coupling      :
*    none.
* --------------------------------------------------------------------------- */

/**
 * @brief Esta función obtiene el mes, dia, hora, minuto y segundo dado el año y los días.
 * @param year Año.
 * @param days días apartir del año year.
 * @param mon Mes. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param day Día. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param hr Hora. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param minute minutos. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param sec segundos. Perderá el valor que tuviera antes de llamar a ésta función.
 */
void    sgp4_date::days2mdhms
        (
          int year, double days,
          int& mon, int& day, int& hr, int& minute, double& sec
        )
   {
     int i, inttemp, dayofyr;
     double    temp;
     int lmonth[] = {31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31};

     dayofyr = (int)floor(days);
     /* ----------------- find month and day of month ---------------- */
     if ( (year % 4) == 0 )
       lmonth[1] = 29;
     i = 1;
     inttemp = 0;
     while ((dayofyr > inttemp + lmonth[i-1]) && (i < 12))
     {
       inttemp = inttemp + lmonth[i-1];
       i++;
     }
     mon = i;
     day = dayofyr - inttemp;

     /* ----------------- find hours minutes and seconds ------------- */
     temp = (days - dayofyr) * 24.0;
     hr   = (int)floor(temp);
     temp = (temp - hr) * 60.0;
     minute  = (int)floor(temp);
     sec  = (temp - minute) * 60.0;
   }  // end days2mdhms

/** -----------------------------------------------------------------------------
*  @brief   this procedure finds the year, month, day, hour, minute and second
*    given the julian date. tu can be ut1, tdt, tdb, etc.
*    @param  jd          - julian date                    days from 4713 bc
*    @param  year       output - year                           1900 .. 2100
*    @param  mon        output- month                          1 .. 12
*    @param  day        output- day                            1 .. 28,29,30,31
*    @param  hr         output- hour                           0 .. 23
*    @param  min        output- minute                         0 .. 59
*    @param  sec        output- second                         0.0 .. 59.999
*    @author  david vallado
*
*
*                           procedure invjday
*
*  this procedure finds the year, month, day, hour, minute and second
*  given the julian date. tu can be ut1, tdt, tdb, etc.
*
*  algorithm     : set up starting values
*                  find leap year - use 1900 because 2000 is a leap year
*                  find the elapsed days through the year in a loop
*                  call routine to find each individual value
*
*  author        : david vallado                  719-573-2600    1 mar 2001
*
*  inputs          description                    range / units
*    jd          - julian date                    days from 4713 bc
*
*  outputs       :
*    year        - year                           1900 .. 2100
*    mon         - month                          1 .. 12
*    day         - day                            1 .. 28,29,30,31
*    hr          - hour                           0 .. 23
*    min         - minute                         0 .. 59
*    sec         - second                         0.0 .. 59.999
*
*  locals        :
*    days        - day of year plus fractional
*                  portion of a day               days
*    tu          - julian centuries from 0 h
*                  jan 0, 1900
*    temp        - temporary double values
*    leapyrs     - number of leap years from 1900
*
*  coupling      :
*    days2mdhms  - finds month, day, hour, minute and second given days and year
*
*  references    :
*    vallado       2007, 208, alg 22, ex 3-13
* --------------------------------------------------------------------------- */

/**
 * @brief Esta función devuelve el año, el mes, el día, la hora, el minuto y el segundo del día juliano pasado por parámetro.
 * @param jd Día juliano.
 * @param year Año. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param mon Mes. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param day Día. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param hr hora. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param minute minuto. Perderá el valor que tuviera antes de llamar a ésta función.
 * @param sec segundo. Perderá el valor que tuviera antes de llamar a ésta función.
 */
void    sgp4_date::invjday
        (
          double jd,
          int& year, int& mon, int& day,
          int& hr, int& minute, double& sec
        )
   {
     int leapyrs;
     double    days, tu, temp;

     /* --------------- find year and days of the year --------------- */
     temp    = jd - 2415019.5;
     tu      = temp / 365.25;
     year    = 1900 + (int)floor(tu);
     leapyrs = (int)floor((year - 1901) * 0.25);

     // optional nudge by 8.64x10-7 sec to get even outputs
     days    = temp - ((year - 1900) * 365.0 + leapyrs) + 0.00000000001;

     /* ------------ check for case of beginning of a year ----------- */
     if (days < 1.0)
       {
         year    = year - 1;
         leapyrs = (int)floor((year - 1901) * 0.25);
         days    = temp - ((year - 1900) * 365.0 + leapyrs);
       }

     /* ----------------- find remaing data  ------------------------- */
     days2mdhms(year, days, mon, day, hr, minute, sec);
     sec = sec - 0.00000086400;
   }  // end invjday


/* -----------------------------------------------------------------------------
*  @brief   this function finds the greenwich sidereal time.
*    @param   jdut1        - julian date in ut1             days from 4713 bc
*    @param   gstime      output- greenwich sidereal time        0 to 2pi rad
*    @author  david vallado
*
*                           function gstime
*
*  this function finds the greenwich sidereal time.
*
*  author        : david vallado                  719-573-2600    1 mar 2001
*
*  inputs          description                    range / units
*    jdut1       - julian date in ut1             days from 4713 bc
*
*  outputs       :
*    gstime      - greenwich sidereal time        0 to 2pi rad
*
*  locals        :
*    temp        - temporary variable for doubles   rad
*    tut1        - julian centuries from the
*                  jan 1, 2000 12 h epoch (ut1)
*
*  coupling      :
*    none
*
*  references    :
*    vallado       2004, 191, eq 3-45
* --------------------------------------------------------------------------- */

/**
 * @brief Esta función devuelve la posición angular del tiempo sidéreo de Greenwich a través de un día juliano.
 * @param jdut1 Día juliano
 * @return el tiempo sidéreo de Greenwich medido en radianes (valores entre 0 t 2pi).
 */
double  sgp4_date::gstime
        (
          double jdut1
        )
   {
     const double twopi = 2.0 * pi;
     const double deg2rad = pi / 180.0;
     double       temp, tut1;

     tut1 = (jdut1 - 2451545.0) / 36525.0;
     temp = -6.2e-6* tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 +
             (876600.0*3600 + 8640184.812866) * tut1 + 67310.54841;  // sec
     temp = fmod(temp * deg2rad / 240.0, twopi); //360/86400 = 1/240, to deg, to rad

     // ------------------------ check quadrants ---------------------
     if (temp < 0.0)
         temp += twopi;

     return temp;
   }  // end gstime


